##faces_analyze.r

##create matrix to store results
tests.matrix = matrix(nrow = 12, ncol = 13)

##function for clustering standard errors
clx <-  function(fm, dfcw, cluster){
  library(sandwich)
  library(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
  u <- apply(estfun(fm),2,
             function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
  coeftest(fm, vcovCL) 
}

##load data
dat = read.csv('faces_data.csv')

##leeop through each test and record results
for(i in 1:5){
  use.dat = dat[dat$test == i,]
  
  ##basic black and white differences from integrated results
  fm.white = summary(lm(white.seg.int~1, data = use.dat))
  fm.black = summary(lm(black.seg.int~1, data = use.dat))
  
  
  ##store these results in matrix
  tests.matrix[i*2,1] = i;   tests.matrix[(i*2)-1,1] = i; 
  
  tests.matrix[(i*2)-1,3] = mean(use.dat$black.int,na.rm=T); tests.matrix[(i*2)-1,5] = mean(use.dat$black.seg,na.rm=T)
  tests.matrix[(i*2)-1,4] = sd(use.dat$black.int,na.rm=T); tests.matrix[(i*2)-1,6] = sd(use.dat$black.seg,na.rm=T)
  
  tests.matrix[i*2,3] = mean(use.dat$white.seg,na.rm=T); tests.matrix[i*2,5] = mean(use.dat$white.int,na.rm=T)
  tests.matrix[i*2,4] = sd(use.dat$white.seg,na.rm=T); tests.matrix[i*2,6] = sd(use.dat$white.int,na.rm=T)
  

  ##extract one-taield p values 
  p.white = pt(coef(fm.white)[, 3], fm.white$df[2], lower=FALSE)  
  p.black = pt(coef(fm.black)[, 3], fm.black$df[2], lower=FALSE)  
  
  tests.matrix[(i*2)-1,7] = fm.black$coefficients[1,'Estimate']; tests.matrix[(i*2)-1,8] = fm.black$coefficients[1,'Std. Error']
  tests.matrix[(i*2),7] = fm.white$coefficients[1,'Estimate']; tests.matrix[(i*2),8] = fm.white$coefficients[1,'Std. Error']

  tests.matrix[(i*2)-1,9] = fm.black$coefficients[1,'t value']
  tests.matrix[i*2,9] = fm.white$coefficients[1,'t value']
  
  tests.matrix[(i*2)-1,10] = fm.black$coefficients[1,'Estimate'] - 1.96*fm.black$coefficients[1,'Std. Error'] 
  tests.matrix[(i*2)-1,11] = fm.black$coefficients[1,'Estimate'] + 1.96*fm.black$coefficients[1,'Std. Error'] 
  
  tests.matrix[(i*2),10] = fm.white$coefficients[1,'Estimate'] - 1.96*fm.white$coefficients[1,'Std. Error'] 
  tests.matrix[(i*2),11] = fm.white$coefficients[1,'Estimate'] + 1.96*fm.white$coefficients[1,'Std. Error'] 
  
  tests.matrix[(i*2)-1,12] = p.black
  tests.matrix[(i*2),12] = p.white
  
  tests.matrix[(i*2)-1,13] = nrow(use.dat)
  
  }

##now pooled results with clustered standard errors by trial
fm.white = lm(white.seg.int~1, data = dat)
fm.white = clx(fm.white,1,dat$test)

fm.black = lm(black.seg.int~1, data = dat)
fm.black = clx(fm.black,1,dat$test)



##store results
tests.matrix[11,3] = mean(dat$black.int,na.rm=T); tests.matrix[11,5] = mean(dat$black.seg,na.rm=T)
tests.matrix[11,4] = sd(dat$black.int,na.rm=T); tests.matrix[11,6] = sd(dat$black.seg,na.rm=T)
tests.matrix[12,3] = mean(dat$white.seg,na.rm=T); tests.matrix[12,5] = mean(dat$white.int,na.rm=T)
tests.matrix[12,4] = sd(dat$white.seg,na.rm=T); tests.matrix[12,6] = sd(dat$white.int,na.rm=T)
tests.matrix[11,7] = fm.black[,'Estimate']; tests.matrix[11,8] = fm.black[,'Std. Error']
tests.matrix[12,7] = fm.white[,'Estimate']; tests.matrix[12,8] = fm.white[,'Std. Error']
tests.matrix[11,9] = fm.black[,'t value']
tests.matrix[12,9] = fm.white[,'t value']
tests.matrix[11,10] = fm.black[,'Estimate'] - 1.96*fm.black[,'Std. Error'] 
tests.matrix[11,11] = fm.black[,'Estimate'] + 1.96*fm.black[,'Std. Error'] 
tests.matrix[12,10] = fm.white[,'Estimate'] - 1.96*fm.white[,'Std. Error'] 
tests.matrix[12,11] = fm.white[,'Estimate'] + 1.96*fm.white[,'Std. Error'] 
tests.matrix[11,12] = fm.black[,'Pr(>|t|)']/2
tests.matrix[12,12] = fm.white[,'Pr(>|t|)']/2

tests.matrix[11,13] = nrow(dat)



##latex versions
tests.matrix = as.data.frame(round(tests.matrix,2))
tests.matrix[,2] = c('Black','White'); tests.matrix[c(11,12),1] = 'Pooled'

tests.matrix[,3] = paste(tests.matrix[,3], ' (',tests.matrix[,4],')',sep='')
tests.matrix[,5] = paste(tests.matrix[,5], ' (',tests.matrix[,6],')',sep='')
tests.matrix[,7] = paste(tests.matrix[,7], ' (',tests.matrix[,8],')',sep='')
tests.matrix[,10] = paste('[',tests.matrix[,10],',',tests.matrix[,11],']',sep='')

tests.matrix[,c(4,6,8,11)] = NULL

colnames(tests.matrix)= c('Trial','Test','Integrated (SD)','Segregated (SD)','Difference (SE)', 'T value',
                               'CI','p','N')


cat('Experiment 1 Results: \n')
print(tests.matrix)

##Print table 1 and s1
table.s1 = xtable(tests.matrix,
                        digits = 2)
table.1 = xtable(tests.matrix[,c(1,2,5:9)],
                  digits = 2)

print(table.1,
      "Table1.tex",
      type = 'latex',
      floating = F,
      booktabs = T)


print(table.s1,
      "TableS1.tex",
      type = 'latex',
      floating = F,
      booktabs = T)


##covariates, table S2
out.matrix = dat[,c('male','non.white','hispanic','us.resident','english',
                                 'college.grad','conservative',
                                 'left.handed','age','income')]

out.matrix.means = as.data.frame(apply(out.matrix,2,mean,na.rm= T))

rownames(out.matrix.means) = c('Male','Non White','Hispanic','US Resident','English Speaker',
                               'College Graduate','Conservative','Left Handed','Age','Income')
colnames(out.matrix.means) = c('Proportion')

table.s2= xtable(out.matrix.means,
                   digits = 2)
print(table.s2,'TableS2.tex',
      type = 'latex',
      floating = F,
      booktabs = T)


